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Abstract 

Energy  conversion  in  electrochemical  power  devices  is  increasingly  correlated  with  high  density  waste  heat  generation,  since 
cell  stacks  have  become  highly  compact.  To  predict  the  transient  thermal  behaviour  of  fuel  cells,  a  three-dimensional  simulation 
of  their  stack  operation  has  been  developed.  Time  constants  of  the  transient  response  following  load  changes  are  determined 
and  are  in  good  agreement  with  experimental  data  obtained  from  a  100  W  fuel  cell  stack. 
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1.  Introduction 

In  recent  years,  the  improvement  of  fuel  cell  power 
density  was  a  development  goal  of  high  priority.  How¬ 
ever,  an  essential  problem  of  high  density  power  devices 
is  the  rejection  of  a  considerable  amount  of  waste  heat 
from  a  small  volume.  Thus,  understanding  the  thermal 
behaviour  of  fuel  cells  during  stationary  and  transient 
operation  receives  increasing  attention.  In  particular, 
knowledge  about  the  temperature  distribution  within 
the  fuel  cell  stack  is  required  to  avoid  thermal  damage 
during  its  operation.  However,  there  is  little  experi¬ 
mental  data  available,  since  probes  cannot  be  integrated 
into  the  compact  electrode  stack.  Thus,  simulation  of 
the  stacks’  thermal  behaviour  is  essential  to  gain  a 
better  understanding  of  heat  generation  and  dissipation 
processes. 

The  present  investigation  is  an  extension  of  our  recent 
work  [1]  on  a  simple  one-dimensional  model  to  predict 
spatial-  and  time-dependent  temperature  profiles  in 
alkaline  fuel  cells  designed  in  the  VARTA  Eloflux 
concept  [2-5], 

The  simulation  presented  in  Ref.  [1]  is  based  on 
rough  simplifications  since  thermal  conductivity  is  re¬ 
stricted  to  the  direction  perpendicular  to  the  electrodes 
and  no  spatial  temperature  distribution  has  been  as¬ 
sumed  across  the  electrodes.  Nevertheless,  these  simple 
calculations  succeed  in  predicting  the  temperature  ver¬ 
sus  load  characteristics  of  an  equilibrium  state  fuel  cell 
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operation  and  give  a  rough  estimate  of  time  constants 
for  thermal  relaxation  after  load  changes.  However, 
the  one-dimensional  model  failed  to  describe  the  de¬ 
tailed  time-dependent  temperature  profiles  during  dy¬ 
namic  fuel  cell  operation. 


2.  Theory 

Considering  a  certain  volume  element  of  the  fuel 
cell,  heat  transfer  across  its  boundaries  and  its  own 
heat  generation  will  cause  changing  temperatures  and 
thus  a  variation  of  the  internal  energy  of  the  volume 
element.  The  time  dependence  of  the  internal  energy 
can  be  expressed  by: 


where  p  and  cp  denote  mass  density  and  the  heat 
capacity,  respectively;  du  describes  the  internal  energy 
of  the  volume  element  unfolded  by  dr,  dy  and  dz. 

The  temperature  of  this  volume  changes  due  to  the 
thermal  flux  across  its  boundaries.  Since  the  heat  con¬ 
ductivity  A  generally  is  assumed  to  be  anisotropic,  the 
partial  derivation  of  the  spatial  temperature  distribution 
is  given  by  three  differential  equations  for  the  three 
dimensions,  respectively: 

a  T 

(2a) 
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vector  f—r0=(xh  ym,  zn)  =  (l Ax,  mAy,  nAz),  and  to  the 
time  tk—t0=k  At,  respectively. 
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where  the  q,  (i=x,  y,  z)  denotes  the  thermal  fluxes. 
Striking  the  thermal  balance  of  this  volume  element, 
including  the  heat  qE  generated  inside,  leads  to  the 
continuity  equation  in  its  differentiated  form: 


+<7e  =  0  (3) 

Considering  those  volume  elements  that  are  run 
through  by  the  coolant  flow,  the  enthalpy  flux  within 
the  coolant  contributes  to  the  thermal  balance.  Thus, 
the  last  term  of  Eq.  (3)  must  be  extended  by  the  thermal 
flux  qF  =  (mFcp,FA  7F)/(dc  dy  dz)  that  is  carried  within 
the  coolant  flow.  Neglecting  the  temperature  depen¬ 
dence  of  the  heat  conductivity  A  simplifies  Eq.  (3)  and 
results  in  an  inhomogeneous,  linear  partial  differential 
equation  of  the  second  order  that  is  known  as  a  Fourier 
differential  equation: 
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Starting  from  a  well-known  temperature  distribution 
as  assumed  for  t=tk  we  can  calculate  the  temperature 
distribution  for  t=tk+1  =4  +  At.  In  a  three-dimensional 
Cartesian  matrix,  unfolded  by  the  parameters  l,  m  and 
n,  the  transient  temperature  distribution  can  be  cal¬ 
culated  by  the  following  approximation: 
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where  the  ah  defined  as  a,  =  A,/cpp  ( i=x,y ,  z),  reflects 
the  anisotropy  of  the  thermal  conductivity. 

Since  the  fuel  cell  is  composed  of  various  layers  of 
different  materials,  an  analytical  solution  of  Eq.  (4) 
cannot  be  given.  Thus,  a  numerical  formulation  of  this 
problem  is  required.  As  an  approach  to  the  exact 
solution,  the  continuous  temperature  distribution  is 
replaced  by  a  discrete  network  of  temperatures. 

A  detailed  description  of  the  difference  method  ap¬ 
plied  to  solve  Eq.  (4)  is  given  in  Ref.  [6]  for  the  one¬ 
dimensional  case.  This  method  is  particularly  suitable 
for  describing  heat-conduction  phenomena  across  mul¬ 
tiple  layers  of  different  materials. 

In  this  discrete  notation  Eq.  (4)  can  be  written  as: 
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In  order  to  determine  a  certain  temperature  T,  m  n  k 
in  a  four-dimensional  parameter  array,  the  temperature 
differences  in  Eq.  (5)  are  expressed  by  the  array  ele¬ 
ments.  Here,  indices  l,  m,  n  and  k  refer  to  the  position 


This  formulation  of  the  heat-conduction  problem 
implies  the  simplification  that  the  balance  of  the  energy 
flows  as  expressed  in  Eq.  (3)  will  not  change  during 
At.  However,  the  thermal  balance  changes  continuously 
and  causes  stability  problems  in  the  numerical  calcu¬ 
lations.  The  finite  spacing  of  the  matrix  elements  (Ax, 
Ay,  Az,  At)  are  crucial  for  the  effort,  precision  and 
stability  of  the  difference  method.  In  particular,  in¬ 
creasing  At  will  cause  thermal  flows  higher  than  needed 
to  even  out  the  temperature  differences  between  neigh¬ 
bouring  volume  elements  and  subsequently  change  the 
sign  of  the  temperature  gradient.  This  will  cause  an 
instability  in  the  difference  method  applied. 

These  oscillations  can  be  avoided  by  introducing  a 
stability  criterion  to  determine  an  upper  limit  for  the 
time  interval  At,  depending  on  the  width  of  the  spatial 
intervals  Ax,  Ay  and  Az.  For  one  dimensional  heat 
conduction  this  stability  criterion  is  defined  assuming 
the  Binder-Schmidt  method  [7]  by:  p=aAt/(Ax)2.  The 
value  of  p  can  be  simply  estimated  assuming  steady- 
temperature  distribution,  and  Tl  k  differs  from  this 
temperature  level  only.  The  restriction  to  one  dimension 
allows  a  temperature  balancing  over  three  neighbouring 
volumes  1—1,1,  and  l+l  only.  Considering  the  limiting 
conditions  of  our  assumption  and  neglecting  heat  sources 
and  sinks  leads  to  a  one-dimensional  expression  of  Eq. 
(7): 
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ATlAk^2pATALk  (8) 

where  AT,tAk  =  Tlk+1-T,tk  and  ATAlk  =  Tt±lk-T/k 
denote  the  temporal  and  spatial  temperature  differ¬ 
ences,  respectively.  Since  the  temporal  temperature 
changes  of  volume  l  must  not  exceed  the  spatial  tem¬ 
perature  differences  to  its  neighbours,  ATl  Ak<ATAUi, 
and  the  stability  criterion  can  be  estimated  to: 

Px  <  1/2  (9) 

in  its  extension  to  three  dimensions  Eq.  (9a)  is  written 
as: 


P*+Py+Pz<  1/2  (9b) 

From  this  we  derive  an  upper  limit  for  At  to  avoid 
oscillations  during  the  numerical  calculation: 


A 


I /_£*-  + -£.  +  -5lV 

2  \Ax2  Ay2  A z2) 


(10) 


Since  this  stability  condition  must  be  fulfilled,  that 
layer  which  leads  to  the  lowest  At  determines  the 
spacing  of  the  grid.  Although  this  model  reflects  the 
real  condition,  it  is  not  practicable  since:  (i)  the  number 
of  volume  elements  is  limited  by  the  capacity  of  the 
computer,  and  (ii)  the  time  required  for  calculation 
increases  by  a  power  of  four  as  the  grid  spacing  de¬ 
creases. 

In  order  to  obtain  moderate  computing  conditions 
with  a  highest  reliability  of  the  results,  heat  conduction 
parallel  to  the  electrode  surfaces  (x-y-plane)  is  treated 
different  from  that  perpendicular  to  the  electrodes  (z- 
direction).  For  the  z-direction,  that  predominates  the 
overall  thermal  balancing,  the  spacing  of  At  is  deter¬ 
mined  according  to  Eq.  (10).  Thermal  balancing  in  the 
x  -y  plane  is  considered  subsequently  to  the  heat  con¬ 
duction  directed  parallel  to  its  normal  vector.  To  con¬ 
tribute  to  tolerable  computing  requirements,  only  a 
rough  spacing  of  Ax  and  Ay  is  considered.  This  coarse 
spacing  within  the  structural  elements  causes  a  re¬ 
alignment  of  the  heat  conductivity  A,  and  ky  and  thus 
a  re-alignment  of  the  factors  ax  and  ay: 


i  -  —p 
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Note  that  this  simplification  will  not  affect  the  heat 
conduction  from  the  source  to  the  sink,  since  they  are 
located  in  different  structural  elements. 

In  the  present  simulation,  we  consider  different  con¬ 
tributions  to  waste-heat  generation,  dissipation  and 
rejection: 

•  reversible  and  irreversible  heat  generation  during  the 
electrochemical  reaction:  entropy  changes;  polari¬ 
zation  losses,  and  ohmic  losses  within  the  electrolyte 


•  Joulean  heat  generation  within  the  electrodes 

•  heat  transition  from  the  structural  elements  to  the 

electrolyte  (coolant)  flow 

•  free  convection  at  the  fuel  cell’s  surfaces 

The  waste  heat  generated  during  fuel  cell  operation 
has  sufficiently  been  discussed  in  the  past  [5,8,9],  Its 
amount  is  mainly  determined  by  the  entropy  changes 
and  load-dependent  overpotentials.  Thus,  the  heat  gen¬ 
erated  can  be  given  by  Eq.  (12)  as  a  function  of  applied 
current  density  and  temperature  at  the  reaction  site: 

0(«rc,  T)  ~  AH—zFufc(itc,  T) 

=  AH— zF[U*  -  r0i(c  &xp(EJRT}]  (12) 

where  Q,  ufc,  i[c  and  T  denote  the  waste-heat  generation 
rate,  the  fuel  cell’s  voltage,  current  density,  and  tem¬ 
perature,  respectively.  The  parameters  U*,  r0  and  Ea 
denote  the  open-loop  voltage,  the  ohmic  resistance  and 
the  activation  energy  of  the  fuel  cell  system  under 
consideration.  These  parameters  depend  on  individual 
cell  configuration  and  geometry  and  can  be  derived 
experimentally. 

The  finite  ohmic  resistance  within  the  conducting 
grid  of  the  electrodes  leads  to  a  potential  drop  across 
the  electrode  and  thus  to  local  potentials  ulc=u{c  (x, 
y)  and  current  densities  itc =i[c  (x,y).  The  spatial  potential 
distribution  can  be  calculated  according  KirchhofFs  law. 
This  potential  distribution  leads  to  inhomogeneous  heat 
generation  and  thus  to  an  inhomogeneous  temperature 
profile  across  the  electrodes’  surfaces.  This  inhomo¬ 
geneity  of  the  temperature  feeds  back  onto  the  waste- 
heat  production  rate  via  the  temperature  dependence 
of  the  fuel  cell’s  characteristic  (see  Eq.  (12)).  Since 
the  ohmic  resistance  is  generally  temperature  depen¬ 
dent,  the  feedback  of  the  temperature  distribution  on 
to  the  local  potential  must  be  considered. 

In  addition,  Joulean  heat  generation  within  the  elec¬ 
trodes  contributes  to  the  inhomogeneity  of  the  tem¬ 
perature  profile  due  to  the  anisotropic  current-collection 
geometry. 

Heat  transition  from  the  structural  materials  to  the 
coolant  flow  and  ambience  has  been  considered  ac¬ 
cordingly  to  the  appropriate  Nusselt  relations: 


Nu=f(Re,Pr)=~d  (13a) 

and 

Nu=f(Gr,Pr)=  (13b) 

for  the  heat  transition  to  the  coolant  drain  and  for 
free  convection,  respectively.  Here  Re,  Pr  and  Gr  denote 
the  Reynolds,  Prandel  and  Grashof  numbers,  respec¬ 
tively.  A  detailed  description  of  the  Nusselt  relations 
was  recently  given  in  Refs.  [1,10]. 
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3.  Experimental 

Experiments  were  carried  out  in  order  to  obtain  the 
characteristic  parameters  of  the  fuel  cell  stack  consid¬ 
ered  as  reflected  in  Eq.  (12).  On  the  other  hand,  the 
simulation  has  to  be  verified  by  experimental  data.  For 
this,  the  fuel  cell  stack  has  been  modified  during 
production  to  enable  the  measurement  of  distinct  tem¬ 
peratures  from  inside  the  stack.  The  fuel  cell  consists 
of  several  layers  performing  different  tasks: 

•  the  electrodes  serve  for  the  electronic  conduction 
and  the  chemical  reactions  taking  place  on  their 
surfaces 

•  the  separators  insulate  the  electrodes  from  each  other 
and  are  filled  with  an  electrolyte  to  enable  ionic 
conduction 

•  the  so-called  ‘labyrinth  layer’  serves  as  a  feed  through 
for  electrolyte  and  coolant  flows 

Due  to  the  compact  arrangement  of  the  stack  layers 
(e.g.,  electrodes  and  separators)  it  is  not  possible  to 
place  temperature  probes  on  the  electrodes  and  si¬ 
multaneously  ensure  safe  fuel  cell  operation.  Thus, 
experimental  temperature  profiles  are  not  obtainable 
from  the  electrodes. 

The  peripheral  experimental  setup  serves  the  coolant/ 
electrolyte  management  and  the  reactant  supply  of  the 
fuel  cell.  The  electrolyte  circuit  consists  of  an  expansion 
tank  containing  the  temperature  control,  a  pump  to 
revolve  the  electrolyte,  and  the  casing.  The  product 
water  is  removed  by  evaporation  and  subsequent  con¬ 
densation  on  cooling  traps  of  a  secondary  coolant  circuit. 
The  expansion  tank  simultaneously  serves  as  evaporating 
tray  for  this  water  removal  process. 

The  reactant-gas  supply  contains  a  humidification 
device  to  avoid  damaging  of  the  fuel  cell  by  drying  out 
their  electrodes.  Inert  gases  that  accumulate  during 
operation  within  the  fuel  cell  can  be  removed  by  purging. 

Considering  further  physical  dimensions  contributing 
to  the  thermal  management  of  the  fuel  cell,  the  ex¬ 
perimental  setup  was  equipped  with  pressure  gauges 
and  flow  meter  for  both  electrolyte  and  reactant  gases. 
In  addition,  the  cell  voltages  and  the  current  generation 
were  monitored.  The  operating  conditions  of  the  fuel 
cell  were  adjusted  by  an  electronic  load. 

The  temperatures  were  recorded  using  NiCr-Ni  ther¬ 
mocouples.  Eight  thermocouples  were  placed  in  each 
labyrinth  layer,  as  indicated  in  Fig.  1.  The  channels 
shown  in  the  Figure  serve  the  electrolyte/coolant  feed. 
The  experiments  were  carried  out  using  a  100  W  fuel 
cell  stack  which  contains  three  labyrinth  layers.  Thus, 
24  temperatures  are  experimentally  available  across  the 
cell  stack. 

In  order  to  determine  the  thermal  response  to  load 
changes  during  fuel  cell  operation,  time-dependent 
temperature  profiles  must  be  recorded.  These  mea¬ 
surements  were  performed  by  a  computer  via  an  AD 


Fig.  1.  Labyrinth  layer  of  the  Eloflux  fuel  cell  with  (A)  coolant  drains 
and  (B)  gas  feedthroughs.  The  positions  of  the  different  thermocouples 
are  indicated  (1  to  8).  All  distances  are  given  in  mm. 

converter  and  an  IEEE  interface.  Since  only  one  ADC 
was  available,  the  data  lines  were  multiplexed.  Thus, 
the  transient  resolution  is  restricted  to  15  s  when 
recording  all  accessible  data. 


4.  Results  and  discussion 

4.1.  Determination  of  the  characteristic  fuel  cell 
parameters 

Consideration  of  the  temperature-dependent  waste- 
heat  generation  in  thermal  modelling  of  cell  stacks 
requires  an  exact  knowledge  of  the  fuel  cell’s  efficiency 
as  a  function  of  temperature  and  current  density.  Thus, 
the  fuel  cell’s  characteristic  was  experimentally  deter¬ 
mined  for  different  temperatures  ranging  from  ambient 
temperature  to  7=60  °C.  These  characteristics  are 
shown  in  Fig.  2;  they  show  almost  linear  behaviour 
except  for  low  current.  From  these  data,  a  Active  open- 
loop  voltage  of  U*  =  1005  mV  is  derived  (see  Fig.  2) 
which  shows  nearly  no  dependence  on  temperature.  In 
accordance  with  Winsel  [5],  the  nearly  linear  section 
of  the  temperature-dependent  fuel  cell  characteristic 
can  be  approximated  by: 

ufc  =  U*  —  r0ifC  cxp(EJRT)  (14) 

where  r0  denotes  a  resistance  characteristic  for  the  fuel 
cell  system  considered  and  Ea  a  characteristic  activation 
energy. 
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Fig.  3.  Experimental  determination  of  fuel  cell  properties  as  activation 
energy  and  internal  resistivity  by  an  Arrhenius’s  evaluation  of  the 
temperature-dependent  characteristics. 


The  coefficients  r0  and  E.d  were  determined  applying 
an  Arrhenius’s  evaluation.  According  to  Eq.  (15),  the 
logarithm  of  the  potential  difference  ln[(t/* -«(/))//]  is 
plotted  against  the  inverse  temperature,  as  shown  in 
Fig.  3;  a  linear  relationship  should  be  found.  Then  the 
activation  energy  Ed  and  the  specific  resistance  r0  can 
be  derived  from  the  slope  and  the  intercept,  respectively: 


-ln(r„)  + 


£.  1 
R  T 


(15) 


For  the  fuel  cell  stack  used  in  our  experiments  we 
found  activation  energies  of  E„  =  20.45  kJ/mol  and 
En  =  21.89  kJ/mol  for  the  two  cell  packages,  respectively. 
The  internal  specific  resistance  could  be  determined 
to  ra  =  85  mfl  cm2  and  r0  =  56  mft  cm2.  These  results 
agree  well  with  those  reported  in  the  literature  for  this 
type  of  fuel  cell  [5], 


4.2.  Equilibrium-state  fuel  cell  temperature 

Since  the  time  constants  of  thermal  relaxation  are 
evaluated  in  terms  of  temperature  differences  with 
respect  to  the  expected  final-state  equilibrium  tem¬ 
perature  (see  Eq.  (19)),  the  characteristic  of  equilibrium 
temperatures  versus  fuel  cell  load  should  be  discussed 
prior  to  the  consideration  of  the  unsteady  case.  Secondly, 
the  inhomogeneity  of  the  spatial  temperature  distri¬ 
bution  will  limit  fuel  cell  operation  at  higher  operating 
temperatures.  We  give  an  attempt  to  evaluate  the  load 
dependence  of  the  inhomogeneity  in  terms  of  fit  pa¬ 
rameters  derived  from  the  experimental  or  calculated 
temperature  versus  load  characteristic. 

The  spatial  distribution  of  equilibrium-state  fuel  cell 
temperatures  has  been  calculated  for  different  operating 
conditions  of  the  cell  stack.  Fig.  4  shows  the  results 
of  this  calculation  for  selected  positions  within  the  cell 
stack.  For  comparability  with  experimental  data,  the 
temperatures  presented  are  related  to  those  volume 
elements  corresponding  to  the  thermocouple  positions 
in  the  experiment.  In  Fig.  4,  the  positions  within  the 
cell  stack  are  indicated  by  the  indices  of  the  corre¬ 
sponding  thermocouples  TC^,.  Here  x  and  y  refer  to 
the  labyrinth  layer  and  to  the  position  within  the  layer 
as  given  in  Fig.  1. 

The  data  are  plotted  against  the  current  generation, 
which  represents  the  fuel  cell’s  load.  Obviously,  the 
equilibrium  temperatures  within  the  fuel  cell  stack  do 
not  depend  linearly  on  the  fuel  cell  load. 

In  order  to  find  a  quantitative  expression  for  fuel 
cell  temperature  as  a  function  of  load,  we  plotted  the 
logarithm  of  AT  versus  the  logarithm  of  the  current 
generation  and  found  a  linear  relationship,  as  also 
shown  in  Fig.  4: 
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(a)  current  generation  in  A  — (b>  current  generation  in  A  — » 

Fig.  4.  Experimentally  derived  equilibrium  temperatures  vs.  fuel  cell 
load.  The  temperatures  shown  were  measured  with  thermocouples 
nos.  2  and  6  in  the  second  cell  of  the  stack,  (a)  The  temperatures 
were  shown  for  different  coolant  flows,  (b)  Comparison  with  calculated 
data.  The  temperature  differences  refer  to  the  equilibrium  temper¬ 
ature  at  /  =  0  A,  i.e.,  the  electrolyte  feed  temperature. 
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ln(AT')  =  ln  a  +  b  Ini' 


(16) 


Here,  AT  denotes  the  difference  between  the  equilib¬ 
rium  temperatures  T(T)  and  T(I=0).  The  latter  is  mainly 
determined  by  the  electrolyte-feed  temperature.  Note 
that  AT'  and  AT  are  nondimensional  quantities  (e.g. 
/'  =///„,  I0=  1  A).  Since  coefficient  b  differs  from  unity, 
the  temperature  does  not  depend  linearly  on  the  current 
generation.  This  behaviour  is  observable  in  both  the 
simulation  results  and  the  experimental  data.  The  ex¬ 
perimentally  derived  values  for  coefficient  b  vary  from 
1.02  1.55.  Note,  that  this  formulation  of  the  tem¬ 

perature  versus  load  characteristic  of  the  fuel  cell  is 
only  a  mathematical  fit  and  is  not  based  on  physical 
principles.  Thus,  parameters  a  and  b  are  not  related 
to  physical  dimensions;  however,  they  may  be  used  to 
compare  the  model  calculation  with  experimental  re¬ 
sults. 

On  the  other  hand,  the  coefficients  from  Eq.  (16) 
can  be  used  to  discuss  the  behaviour  of  the  temperature 
distribution  as  a  function  of  fuel  cell  load.  Increasing 
the  current  density,  diverging  temperatures  at  random 
but  different  positions  within  the  friel  cell  stack  should 
result  in  a  more  inhomogeneous  temperature  distri¬ 
bution.  Since  the  convergence  behaviour  of  the  tem¬ 
perature  distribution  can  be  expressed  by  the  ratio  of 
two  temperatures  rather  than  by  their  difference,  the 
inhomogeneity  of  the  temperature  distribution  can  be 
expressed  by: 


/conv  =  ^ 
J  AT2' 


(17) 


where  /conv,  denotes  the  convergence  factor  and  indices 
1  and  2  refer  to  different  volume  elements  of  the  fuel 
cell  stack.  Convergence  factors  /conv  <  1  will  lead  to  a 
more  homogeneous  temperature  distribution,  increasing 
the  fuel  cell  load,  whereas  /conv  >  1  will  increase  the 
inhomogeneity  of  the  temperature  distribution. 

From  Eq.  (17),  the  temperature  ratios  are  obviously 
not  affected  by  load  variations  for  identical  coefficients 
b,=bj  of  the  data  field.  However,  the  experimental 
results  indicate  a  variation  of  £>,  depending  on  the 
position  of  the  related  volume  element  within  the  stack. 
Thus,  we  conclude,  that  the  temperature  ratios  of  the 
different  volume  elements  i,  j  vary  with  increasing  load. 
Depending  on  the  coefficients  a,  b  of  the  polynomial 
in  Eq.  (16),  the  temperature  ratios  and  thus  the  in¬ 
homogeneity  of  the  temperature  distribution  may  in¬ 
crease  or  decrease  with  increasing  load.  However,  if 
the  fuel  cell  load  is  sufficiently  high,  the  temperature 
ratios  will  diverge,  since  the  polynomial  with  the  higher 
exponent  grows  faster.  These  load  limits  can  be  esti¬ 
mated  by  determining  the  intersections  of  two  random 
polynomials.  Then,  the  limiting  condition,  that  the  fuel 
cell  load  must  exceed  a  lower  limit  to  show  clearly  a 
diverging  temperature  distribution,  can  be  expressed 


for  two  sets  of  parameters  a,  b  by  a  divergence  criterion 
Cd  defined  as: 

/  s  (U'v-ti)) 

r>cH«J  (18) 

Evaluating  coefficients  a,  and  bh  we  found  a  wide 
variation  of  the  defined  divergence  criterion.  The  fre¬ 
quency  distribution  of  this  divergence  criterion  n(Cd) 
is  shown  in  Fig.  5  by  the  dashed  line.  The  integrated 
frequency  distribution  Ar(Cd)  =  /^d„/j(cd)  dcd  as  given 
by  the  solid  line  represents  the  probability  that  the 
divergence  criterion  (Eq.  (18))  is  fulfilled  for  random 
pairs  of  parameter  sets. 

Evaluating  the  statistical  distribution  of  this  diver¬ 
gence  criterion  we  found  that  75%  of  the  possible 
parameter  combinations  fulfil  this  criterion  in  the  load 
range  of  interest  (20  A</<60  A).  Thus  we  conclude 
that  spatial  temperature  distribution  will  diverge  with 
increasing  fuel  cell  load. 

Since  the  I  versus  U  characteristic  of  the  fuel  cells 
depends  on  temperature,  increasing  inhomogeneity  in 
the  thermal  distribution  causes  increasing  inhomo¬ 
geneity  of  the  current  density. 

Although  the  equilibrium  temperatures  do  not-depend 
linearly  on  the  load,  the  curvature  of  the  Teq  versus 
I  characteristic  is  small  in  the  considered  load  range 
(20  A</<50  A)  and  only  four  data  points  are  available 
from  both  dynamic  simulation  calculation  and  exper¬ 
iment.  Due  to  the  experimental  error  bars  and  the  few 
data  points,  the  evaluation  of  the  coefficients  from  Eq. 
(16)  is  difficult  and  not  suitable  for  comparing  the 
system  simulation  with  the  experimental  data. 

However,  as  can  be  seen  from  Fig.  6,  the  Teq  versus 
I  characteristic  can  almost  be  approximated  by  straight 
lines  whose  slope  is  related  to  the  mean  slope  of  the 
nonlinear  expression  of  Eq.  (16)  by  the  relation 


log((ai/aj)1/<brbi))  - <~ 

Fig.  5.  Frequency  distribution  of  the  divergence  criterion  as  defined 
in  Eq.  (18)  (dashed  line)  and  the  integrated  frequency  distribution 
(solid  line),  which  is  interpreted  as  the  probability  of  diverging 
temperature  behaviour  in  the  fuel  cell  stack. 
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A77A/=ai/(f’_1>.  Here  7  denotes  that  current  whose 
derivative  equals  the  mean  slope  of  the  monotone 
function  in  a  given  range. 

Fig.  7  gives  a  comparison  of  the  crude  linear  ap¬ 
proximation  of  the  simulation  data  to  the  more  detailed 
evaluation  using  Eq.  (16).  Assuming  7=35  A,  we  found 
a  good  correlation  between  both  evaluations  as  indicated 
by  the  parameters  m  =  1.037  and  a=-3.9XlO-3  K 
A-1  of  the  regression  line.  Thus,  the  mean  slopes  of 
the  Teq  versus  7  characteristics  are  suitable  for  comparing 
the  experiment  and  the  system  simulation.  This  com¬ 
parison  is  shown  in  Fig.  8.  Although  the  plotted  data 
show  a  distinct  scattering,  the  experimental  data  agree 
sufficiently  with  the  system  simulation  results,  as  can 
be  seen  from  the  regression  line  parameters  m  =  0.96 
and  a  =  0.12  K  A-1. 

From  this  we  conclude  that  the  presented  simulation 
program  is  suitable  for  giving  a  correct  description  of 
the  dependence  of  equilibrium  temperatures  on  fuel 
cell  load. 

4.3.  Thermal  response  on  load  changes 

The  reason  for  performing  dynamic  calculations  of 
heat  and  mass  transport  in  fuel  cells  is  to  gain  infor- 


ai*bi*T<br1>  - ► 

Fig.  7.  Comparison  of  the  crude  linear  evaluation  of  the  Tcq  vs.  I 
characteristic  from  the  system  calculation,  with  the  mean  slope 
determined  from  coefficients  a,,  h,  which  were  obtained  from  a  more 
detailed  nonlinear  evaluation  of  calculated  equilibrium  temperatures. 

mation  on  basic  principles  that  drive  the  thermal  re¬ 
laxation  of  a  fuel  cell  system  after  load  changes.  In 
addition,  we  given  an  attempt  to  relate  the  basic  prin¬ 
ciples  to  the  experimentally  observable  quantities.  The 
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(AT/AI)calc  - ► 


Fig.  8.  Comparison  of  the  experimental  and  calculated  slopes  of  the 
Tcq  vs.  I  characteristics. 


operating  time  in  s  - *- 

Fig.  9.  Thermal  relaxation  of  the  fuel  cell  stack  after  load  changes. 
The  experimental  data  were  measured  with  thermocouples  nos.  2 
and  6  in  the  second  labyrinth  layer.  In  addition,  the  temperature 
of  the  coolant  (KOH)  is  shown. 

basic  principles  of  thermal  relaxation  are  in  general 
the  heat  fluxes  and  the  heat  transition,  as  discussed 
above  in  Section  2.  Experimental  observables  of  those 
processes  are  their  time  constants.  Further,  it  is  im¬ 
portant  to  distinguish  different  contributions  to  the 
relaxation  process  and  to  identify  the  basic  principles 
by  which  they  are  driven. 

The  thermal  response  of  the  fuel  cell  on  load  changes 
was  measured  starting  from  equilibrium-state  fuel  cell 
operation  at  /= 5  A,  corresponding  to  i  =  12.5  mA  cm-2 
and  an  electric  power  generation  of  .Pei  =  10  W.  Tem¬ 
perature  data  were  recorded  from  the  thermocouples 
in  the  labyrinth  layer  after  increasing  the  electric  load 
to  Fei  =  90  W  corresponding  to  7=60  A  or  i  =  150  mA 
cm-2.  Initially  a  steep  increase  of  the  temperature  is 
observed,  as  shown  in  Fig.  9;  however,  the  temperatures 
show  asymptotic  behaviour  a  few  seconds  later.  Equi¬ 
librium  temperatures  have  almost  been  achieved  after 
a  few  hundred  seconds.  After  the  fuel  cell  has  been 


operated  for  1600  s  at  constant-current  generation  of 
7=60  A,  the  load  was  decreased  to  7=20  A,  corre¬ 
sponding  to  i'  =  50  mA  cm-2  or  7^,  =  37.5  W. 

On  the  other  hand,  dynamic  temperature  calculations 
were  performed  for  those  volume  elements  that  cor¬ 
respond  with  the  positions  of  the  thermocouples.  The 
start  parameters  of  this  dynamic  system  modelling  cor¬ 
respond  with  the  initial  conditions  of  the  experimental 
setup.  Running  down  the  calculation,  the  load  parameter 
was  changed  in  order  to  simulate  the  experimental  load 
profile  of  the  fuel  cell. 

In  this  way,  we  generate  two  sets  of  transient  tem¬ 
perature  data,  by  experiment  and  system  simulation. 
These  data  sets  correspond  to  equal  operating  conditions 
of  the  fuel  cell  stack  and  are  suitable  for  validating 
the  system  calculation. 

In  order  to  derive  an  analytical  expression  for  the 
time  dependence  of  the  fuel  cell  temperatures,  we 
suggest  that  the  relaxation  process  after  load  changes 
can  be  qualitatively  described  in  good  approximation 
by  an  exponential  behaviour.  Then  the  transient  coef¬ 
ficients  describing  the  relaxation  process  can  be  eval¬ 
uated  by  plotting  temperature  differences  to  the  final- 
state  equilibrium  temperature  A  T=Tcq  {-T(t)  in  a 
logarithm  scale  versus  time,  as  shown  in  Fig.  10. 

It  is  clearly  seen  that  the  simple  exponential  de¬ 
scription  of  the  transient  temperature  behaviour  failed, 
since  the  plot  did  not  result  in  a  linear  relation  of 
In  AT  and  t.  However,  from  Fig.  10  we  conclude  that 
the  thermal  relaxation  can  be  described  in  general  by 
a  sum  of  different  exponential  terms: 

T(t)  =  Teq,  f-  2  a,  exp(  -  ^t)  (19) 

where  n  denotes  the  number  of  terms,  r,  denotes  the 
time  constants,  and  the  total  difference  between  the 
initial-  and  final-state  equilibrium  temperatures  is  given 
by  Xa/=Teq  f— req>/.  Evaluating  the  data  shows  that 
only  two  terms  («  =  2)  lead  to  a  sufficiently  good  de- 
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scription  of  the  thermal  response  on  load  changes.  From 
this  we  conclude  that  two  processes  contribute  to  the 
thermal  relaxation  of  fuel  cells.  The  fast  process  shows 
transient  constants  in  the  range  of  4.1  XlO-3 
s_1  <r,<  1.8x  10  2  s“\  whereas  the  slower  process 
shows  transient  constants  ranging  from  t,<8.8x10~4 
s_1  to  t,<6.8x1CT3  s_1. 

The  suggestion  that  the  fast  process  may  be  related 
to  heat  conduction  perpendicular  to  the  electrode  sur¬ 
faces  (along  the  z-axis  in  the  simulation  model)  and 
the  slower  process  assigned  to  heat  conduction  within 
the  structural  components  (e.g.,  electrode,  separators), 
which  corresponds  to  the  x-y-plane,  is  based  on  the 
anisotropy  of  the  heat  conductivity.  On  the  other  hand, 
the  relaxation  process  is  very  complex.  Temperature 
changes  affect  the  heat  generation  and  cause  a  re¬ 
alignment  of  the  current-density  distribution.  The  slower 
process  may  reflect  a  re-alignment  process  occurring 
subsequently  to  the  initial  temperature  changes.  How¬ 
ever,  in  our  opinion,  re-alignment  processes  run  down 
quickly  compared  with  the  heat  conduction  within  the 
cell  stack.  Secondly,  the  suggested  anisotropic  behaviour 
of  the  heat  conduction  is  supported  by  the  evidence 
of  only  one  contribution  to  the  thermal  relaxation  from 
a  simple  one-dimensional  thermal  simulation  [1].  Thus 
the  simple  one-dimensional  description  of  the  transient 
temperature  behaviour  failed. 

Temperature  changes  in  a  certain  volume  element 
are  caused  by  heat  fluxes  across  its  boundary,  as  long 
as  no  heat  sources  or  sinks  are  located  in  the  volume 
element  considered.  Changes  in  the  heat  content  of 
the  volume  elements  can  be  expressed  by  the  derivative 
of  Eq.  (19): 

=  -  pcp  dV^oti  Tj  exp(  -  tJ)  =  2(2/  (20) 

On  the  other  hand,  these  heat  fluxes  are  driven  by 
temperature  gradients.  Considering  anisotropic  heat 
conductivity  (i.e.,  A  =  \(<p,i)))  leads  to  the  following 
expression  of  Eq.  (4): 

?  =  / / x(<p’  d)  aI*aS*°  d<p  (21) 

Finally,  from  a  comparison  of  Eqs.  (20)  and  (21) 
and  substitution  of  dT/dt  in  Eq.  (21)  by  the  appropriate 
expression  from  Eq.  (20),  we  conclude  that  the  ex¬ 
perimentally  identified  contributions  to  the  thermal 
relaxation  of  the  fuel  cell  stack  are  related  to  anisotropic 
thermal  conduction  within  the  cell  stack.  Since  we  found 
evidence  of  only  two  contributions  to  the  thermal 
relaxation,  the  anisotropy  of  the  heat  conduction  is 
restricted  to  two  nonequivalent  directions.  To  account 


exp.  time  constants  (s'1)  - <*- 

Fig.  11.  Comparison  of  calculated  and  experimentally  derived  time 
constants  describing  the  thermal  relaxation  of  fuel  cells  during  dynamic 
operation. 


for  two  linear  independent,  nonequivalent  directions, 
Eq.  (21)  should  be  expressed  using  cylindrical  coor¬ 
dinates  (z,  r,  <p).  Then  the  empirically  derived  expression 
in  Eq.  (19)  can  be  related  to  the  heat  conductivity  by 
the  following  relation: 


2>,t,  exp(-,,)--  ^ 


2K#T 
pcp  dr2 


(22) 


To  validate  the  transient  behaviour  of  the  simulation 
model,  the  time  constants  from  Eq.  (19)  have  been 
evaluated  for  both  calculated  and  experimental  data. 
In  Fig.  11,  the  calculated  time  constants  are  plotted 
against  the  experimental  data.  From  a  linear  fit  to  these 
data,  we  derived  a  slope  of  m  =  1.012  and  an  intercept 
of  a=5.1xl0-4  s_1  for  the  regression  line.  Thus,  a 
comparison  of  experimental  and  calculated  time  con¬ 
stants  results  in  a  functional  relation  that  is  very  close 
to  the  identity  operator.  From  this  we  conclude  that 
the  presented  simulation  model  is  suitable  for  giving 
a  correct  description  of  the  time  dependence  of  spatial 
temperature  distributions  with  fuel  cell  stacks  during 
dynamic  operation. 


5.  Conclusions 

A  three-dimensional  time-dependent  simulation  pro¬ 
gram  has  been  developed  in  order  to  describe  the 
temperature  distribution  within  a  fuel  cell  stack  and 
its  response  to  load  changes. 

The  results  of  the  model  calculations  have  been 
verified  by  experimentally  accessible  temperatures  ob¬ 
tained  from  measurements  on  an  Eloflux  fuel  cell  block. 
From  the  comparison  of  experimental  and  calculated 
temperatures  we  found  a  good  agreement  of  the  thermal 
modelling  with  the  experiment  concerning  the  differ¬ 
ences  of  equilibrium  temperatures  related  to  different 
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operating  states  of  the  fuel  cell.  From  the  temperature 
versus  load  characteristic,  we  predict  an  increasing 
inhomogeneity  of  temperature  distribution  with  in¬ 
creasing  fuel  cell  load. 

In  addition,  the  transient  behaviour  of  the  experi¬ 
mental  temperatures  is  described  correctly  by  the  sim¬ 
ulation.  We  found  two  contributions  to  the  thermal 
relaxation  following  the  load  changes  during  stack  op¬ 
eration: 

(i)  a  fast  process  showing  transient  constants  4.1  X 
10-3  s-1  <r<  1.8X  10-2  s-1,  and 

(ii)  a  slower  process  with  transient  constants  in  the 
range  8.8X10-4  s-1<t<6.8x10-3  s-1. 

These  thermal  relaxation  processes  are  related  to 
the  heat  conduction  perpendicular  to  the  electrodes’ 
surfaces  and  to  in-plane  heat  conduction  of  the  structural 
elements.  Although  the  transient  constants  vary  over 
more  than  one  order  of  magnitude,  we  found  a  good 
correlation  between  experimental  and  calculated  tem¬ 
peratures.  Thus,  the  transient  thermal  behaviour  of  the 
fuel  cell  stack  is  correctly  described  by  the  simulation 
program. 


6.  List  of  symbols 

ax ,  ytZ  thermal  conductivity  in  x-,  y-  and  z-directions 

(m2  s-1) 

cd  divergence  criterion 

cp  specific  heat  capacity  (kJ  kg-1  K-1) 

£a  activation  energy  (kJ  mol-1) 

F  Faraday’s  constant  (A  s  mol-1) 

/conv  convergence  factor 

Gr  Grashof  number  (-) 

H  enthalpy  (J  mol-1) 

/  current  (A) 

itc  current  density  of  the  fuel  cell  (mA  cm-2) 

mF  mass  flow  of  the  coolant  (kg  s-1) 

Nu  Nusselt  number  (-) 

P,  Px.y.z  stability  criterion  (-) 

Pr  Prandtl  number  (-) 

qB  heat  generation  rate  (W  m-2) 

qF  thermal  flow  within  the  coolant  (W  m-2) 


qx,y,x  thermal  flow  in  x-,  y-  and  z-directions  (W 
m-2) 

R  universal  gas  constant  (J  mol  1  K  ') 

r0  internal  area  specific  source  resistance  of 

the  fuel  cell  (fl  cm2) 

Re  Reynolds  number  (  — ) 

t  time  (s) 

T  temperature  (K) 

Teq  equilibrium  temperature  (K) 

T,mnk  temperature  of  the  volume  element  /,  m,  n 

at  the  time  k  (K) 
u  internal  energy  (J) 

u[c  voltage  of  the  fuel  cell  (V) 

U*  Active  open-loop  voltage  (V) 

V  volume  (m3) 

Greek  letters 

a  heat- transition  coefficient  (W  m-2  K-1) 

A XtyiX  heat  conductivity  x-,  y-  and  z-directions 

p  density  (kg  m-3) 

t,  time  constant  (s-1) 
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